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A high-significance detection of non-Gaussianity in the WMAP 
1-year data using directional spherical wavelets 



o 
o 

(N 

en 
> 

o 

o 
o 

o 



^ 



J. D. McEwen,^* M. P. Hobson,^ A. N. Lasenby^ and D. J. Mortlock^ 

^Astrophysics Group, Cavendish Laboratory, Madingley Road, Cambridge CBS OHE, UK 
Institute of Astronomy, Madingley Road, Cambridge CBS OHA, UK 



1 February 2004 



ABSTRACT 

A directional spherical wavelet analysis is performed to examine the Gaussianity of the 
Wilkinson Microwave Anisotropy Probe (WMAP) 1-year data. Such an analysis is facilitated 
by the introduction of a fast directional continuous spherical wavelet transform. The di- 
rectional nature of the analysis allows one to probe orientated structure in the data. Sig- 
nificant deviations from Gaussianity are detected in the skewness and kurtosis of spher- 
ical elliptical Mexican hat and real Morlet wavelet coefficients for both the WMAP and 
ISgmarketal. (2003) foreground-removed maps. The previous non-Gaussianity detection 
made by Vielva et al. ( 2003) using the spherical symmetric Mexican hat wavelet is confirmed, 
although their detection at the 99.9% significance level is only made at the 95.3% significance 
level using our most conservative statistical test. Furthermore, deviations from Gaussianity in 
the skewness of spherical real Morlet wavelet coefficients on a wavelet scale of 550' (cor- 
responding to an effective global size on the sky of ~ 26° and an internal size of ~ 3°) at 
an azimuthal orientation of 72°, are made at the 98.3% significance level, using the same 
conservative method. The wavelet analysis inherently allows one to localise on the sky those 
regions that introduce skewness and those that introduce kurtosis. Preliminary noise analysis 
indicates that these detected deviation regions are not atypical and have average noise disper- 
sion. Further analysis is required to ascertain whether these detected regions correspond to 
secondary or instrumental effects, or whether in fact the non-Gaussianity detected is due to 
intrinsic primordial ffuctuations in the cosmic microwave background. 

Key words: cosmic microwave background - methods: data analysis - methods: numerical 



1 INTRODUCTION 

A range of primordial processes may imprint signatures on the 
temperature fluctuations of the cosmic microwave background 
(CMB). The currently favoured cosmological model is based on 
the assumption of initial fluctuations generated by inflation. In the 
simplest inflationary models, these result in Gaussian tempera- 
ture anisotropics in the CMB. Non-standard inflationary models 
and various cosmic defect scenarios could, however, lead to non- 
Gaussian primordial CMB fluctuations. Non-Gaussianity may also 
be introduced by secondary eff'ects, such as the reionisation of the 
Universe, the integrated Sachs-Wolfe eff'ect, the Rees-Sciama ef- 
fect, the Sunyaev-Zel'dovich eff'ect and gravitational lensing - in 
addition to measurement systematics or foreground contamination. 
Consequently, probing the microwave sky for non-Gaussianity is 
of considerable interest, providing evidence for competing scenar- 
ios of the early Universe and also highlighting important secondary 
sources of non-Gaussianity and systematics. 

Ideally, one would like to localise any detected non-Gaussian 
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components on the sky, in particular to determine if they correspond 
to secondary eff'ects or systematics. The ability to probe diff'erent 
scales is also important to ensure non-Gaussian sources present 
only on certain scales are not concealed by the predominant Gaus- 
sianity of other scales. Wavelet techniques are thus a perfect can- 
didate for CMB non-Gaussianity analysis, since they provide both 
scale and spatial localisation. In addition, directional wavelets may 
provide further information on orientated structure in the CMB. 

Wavelets have already been us ed to analyse the G aus- 
sianity of the CMB. Fo r example, iHobson et alJ il999h and 
iBarreiro & Hobsoa i200ll) investigated the use of planar wavelets 
in detecting and characterising non-Ga ussianity on patc hes of the 
CMB sky. This approach was used bv lMukheriee et alJ (2000) to 
analyse planar faces of the 4-year Cosmic Background Explorer- 
Diff'erential Microwave Radiometer (COBE-DMR) data in the 
QuadCube pixelisation, showing that the data is consistent with 
Gaussianity (corre cting an earlier claim of non-Gaussianity by 
IPando et all Il998h . To consider a full sky CMB map properly, 
however, wavelet analysis must be extended to spherical ge- 
ometry. A spherical H aar wavelet analysis of the COBE-DMR 
data was performed by IBarreiro et alJ i2000h . but no evidence of 
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no n-Gaussianity w as found. Employ ing the approach described 
by lAntoine & Vand erghevnsfl il998|) for pe rforming continuous 



wavelet transforms on a sphere, ICavon et alJ (2001) used the 
isotropic Mexican hat wavelet to analyse the COBE-DMR maps; 
again, no significant dev i ations from Gaussianity were detected. 
iMartinez-Gonzalez et alJ i200A subsequently compared the per- 
formance of spherical Haar and Mexican hat wavelets for non- 
Gaussianity detection and found the Mexican hat wavelet to be su- 
perior. 

Since the release of the Wilkinson Microwave Anisotropy 
Probe (WMAP) 1-year data, a wide range of Gaussianity anal- 
yses have been performed, calculating measures such as the 
bispectrum and M inkowski function al ( Komatsu et alJ l2003l : 
iMagueiio & Medeiros 2004.: L and & Magueiio..2004) . the genus 
JCoUev & GotJl2003lilEriksen et alJl200A . the 3-Doint correlation 



function ("Gaztanaga & Wagg"2003'), multipole alignment statistics 
(CoDietal. 2004; de Oliveira-Costa et al. 2004; Slosar & Seliak 
|2004|), phase associations ( Chiang etal. 2003; Coles et al. i2004l) . 
local curvature ( H ansen et al. ■2004 jj^Cabella_etal.. .2004 ) and hot 
and cold spot statistics JLarson & Wandelt|l2004h . Some statistics 
show consistency with Gaussianity, whereas others provide some 
evidence for a non-Gaussian signal and/or an asymmetry between 
the northern and southern Galactic hemispheres. One of the highest 
si gnificance levels for non-Gaussianity yet reported was obtained 
bv lVielva et alJ yOOj) using a spherical Mexican hat wavelet anal - 
ysis. This result has been confirmed by Mukherie e & Wang|J2004i) . 
who show it to be robust to diff'erent Galactic masks and assump- 
tions regarding noise properties. In particular, it was found that the 
kurtosis of the wavelet coefficients in the southern hemisphere, at 
an approximate size on the sky of 10°, lies just outside the 3cr Gaus- 
sian confidence level. 

Previous wavelet analyses of the CMB have been restricted to 
rotationally symmetric wavelets. A directional analysis on the full 
sky has previously been prohibited by the computational infeasi- 
bility of any implementation. In this paper, by applying a fast di- 
rectional continuous spherical wavelet transform (CSWT), we ex- 
tend non-Gaussianity analysis to examine directional structure in 
the CMB. 

The remainder of this paper is structured as follows. The di- 
rectional CSWT and the construction of new directional spherical 
wavelets is presented in Section |5| In Section I^Jthe procedure fol- 
lowed to analyse the WMAP 1-year data for non-Gaussianity is 
described. Results and further analysis are presented in Section |3] 
Concluding remarks are made in Section [S] 



2 DIRECTIONAL CONTINUOUS SPHERICAL 
WAVELET ANALYSIS 

To perform a wavelet analysis of full sky maps defined on the celes- 
tial sphere, Euclidean wavelet analysis must be extended to spher- 
ical geometry. We consider the d irectional CSWT constructed by 
lAntoine & VanderghevnstI il998h . This transform was constructed 
from group theoretic principles, however we present here an equiv- 
alent construction based on a few simple operations and norm- 
preserving properties. 



2.1 Transform 

A wavelet basis is constructed on the sphere by applying the 
spherical extension of Euclidean motions and dilations to mother 




Figure 1. Stereographic projection of the sphere onto the plane. 

wavelets defined on the sphere - analogous to the construction of a 
Euclidean wavelet basis. 

The natural extension of Euclidean motions on the sphere are 
rotations. These are characterised by the elements of the rotation 
group S0(3), which we parameterise in terms of the three Euler 
angles (a,/3,y). The rotation of a square-integrable function / on 
the 2-sphere S^ (i.e. / e L^(S^)) is defined by 



{n,f){aj)^f(p-'aj), peS0(3), 



(1) 



where oj denotes spherical coordinates (i.e. (x) g S^). 

Dilations on the sphere are constructed by first lifting 5 ^ to the 
plane by a stereographic projection from the south pole (Figure [l), 
followed by the usual Euclidean dilation in the plane, before re- 
projecting back onto S^. A spherical dilation is thus defined by 



(DafXco) = faioj) = ^/A(^ fioj.la), «' 



(2) 



where oja = (Oa,(p) and tan(^a/2) = atan(^/2). The A(a,6) cocycle 
term is introduced to preserve the 2-norm and is defined by 



^('^'^^ [(a2-l)cos6' + (a2 + l)]2 ' 

A wavelet basis on the sphere may now be constructed by ro- 
tations and dilations of an admissible^ mother spherical wavelet 
i// e L^(S^) (described further in Section IT2l . The correspond- 
ing wavelet family {if/a,p = KpDa^, p e S0(3), a e R+} provides an 
over-complete set of functions in L^{S^). The CSWT of 5 e L^{S^) 
is given by the projection onto each wavelet basis function in the 
usual manner. 



(3) 



where the ^ denotes complex conjugation and di\i{(jS) - sin(^) d^d0 
is the usual rotationally invariant measure on the sphere. 

The transform is general in the sense that all orientations in 
the rotation group S0(3) are considered, thus directional structure 
is naturally incorporated. It is important to note, however, that only 
local directions make any sense on 5^. There is no global way of 
defining directions on the sphere^ - there will always be some sin- 
gular point where the definition fails. 

A full directional wavelet analysis on the sphere has previ- 
ously been prohibited by the computational infeasibility of any im- 
plementation. We rectify this problem by presenting a fast algo- 
rithm in Appendix lAlto perform the directional CSWT. 

^ Candidate mother wavelets must a lso satisfy certain adm issibility criteria 
to qualify as a spherical wavelet (see lAntoine et al.l200a for a definition of 
the strict admissibility criterion and the more practical, necessary but not 
sufficient, zero mean criterion). 

2 There is no differentiable vector field of constant norm on the sphere and 
hence no global way of defining directions. 
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2.2 Mother spherical wavelets 

The wavelet basis previously described is constructed from ro- 
tations and dilations of an admissible mother spherical wavelet. 
Mother spherical wavelets are simply constructed by projecting ad- 
missible Euclidean planar wavelets onto the sphere by an inverse 
stereographic projection, 

2 



jA,2(^,0) = (n->K2)(^,0) = 



1 + cos(6') 



il/R2(r,(p), 



(4) 



where r = 2tan(^/2). The modulating term is again introduced to 
preserve the 2-norm. 

Directional spherical wavelets may be naturally constructed in 
this setting - they are simply the projection of directional Euclidean 
planar wavelets onto the sphere. Two directional planar Euclidean 
mother wavelets are defined in the following subsections: the el- 
liptical Mexican hat and real Morlet wavelets. The corresponding 
spherical wavelets are illustrated in Figure |2| The Mexican hat 
wavelet and the real Morlet wavelet, chosen for its sensitivity to 
scanning artifacts, are subsequently applied to the detection of non- 
Gaussianity in the WMAP 1-year data. 



2. 2. 1 Elliptical Mexican hat wavelet 

We propose a directional extension of the usual Mexican hat 
wavelet. The elliptical Mexican hat wavelet is defined as the nega- 
tive of the Laplacian of an elliptical 2-dimensional Gaussian, 
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which reduces to the usual symmetric Mexican hat wavelet for the 
special case where cr^ - cTy. The elliptical Mexican hat wavelet is 
invariant under integer azimuthal rotations of n, thus the rotation 
angle y is always quoted in the range [0, n). 

We define the eccentricity of an elliptical Mexican hat wavelet 
as the eccentricity of the ellipse defined by the first zero-crossing, 
given by 



1-1^ 



(6) 



Elliptical Mexican hat wavelets are subsequently parameterised by 
their eccentricity; the standard deviation in each direction is set by 



(Tv = 1 and (Tv 



,vr 



6^. Elliptical Mexican hat wavelets are 



illustrated in Figure |2| (a) and (b) for eccentricities e - 0.00 and 
6 = 0.95 respectively. 

We define the eff'ective size on the sky of a spherical ellip- 
tical Mexican hat wavelet for a particular dilation as the angular 
separation between the first zero-crossings on the major axis of the 
ellipse, given by 



^r(«) 



■ 4 tan 



4U2V2a. 

V2/ 



(7) 



2. 2. 2 Real Morlet wavelet 

The real Morlet wavelet is defined by 



il^\x;k) = Kt 



exp 



ik X 

~1T 



exp 



(8) 



where k is the wave vector of the wavelet (henceforth we consider 
only wave vectors of the form k - (ko,0)^). We have scaled the 
usual definition of the real Morlet wavelet to achieve size consis- 
tency with the elliptical Mexican hat wavelet. The real Morlet 
spherical wavelet is also invariant under integer azimuthal rotations 
of ;7r, thus the rotation angle y is always quoted in the domain [0, tt). 
The real Morlet wavelet has two orthogonal scales: one defin- 
ing the overall size of the wavelet and the other defining the size 
of its internal structure. The overall eff'ective size on the sky of the 
spherical real Morlet wavelet is defined as the angular separation 
between opposite e~^ roll-off" points of the exponential decay fac- 
tor, and is given by 



^f°^(a) = 4tan- 



a 



2V2a. 



(9) 



Notice that for a given dilation a, the spherical elliptical Mexican 
hat and real Morlet wavelets have an equivalent overall eff'ective 
size on the sky. The eff'ective size on the sky of the internal struc- 
ture of the real Morlet wavelet is defined as the angular separation 
between the first zero-crossings in the direction of the wave vector 
k, and is given by 



^T\^) ^ 



4 tan 
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ko 



(10) 



3 NON-GAUSSIANITY ANALYSIS 

Spherical wavelet analysis is applied to probe the WMAP 1-year 
data for possible deviations from Gaussianity. We follow a similar 
strategy to Vielva et al. (2003), however we extend the analysis to 
directional spherical wavelets to probe orientated structure in the 
CMB. 



3.1 Data preprocessing 

We co nside r the same data se t analysed by both iKomatsu et all 
i2003h and I Vielva et all i2003h in their non- Gaussianity studies. 
The observed WMAP maps for which the CMB is the domi- 
nant signal (two Q-band maps at 40.7GHz, two V-band maps at 
60.8GHz and four W-band maps at 93.5GHz) are combined to 
give a single signal-to-noise ratio enhanced map. These maps, to- 
gether with receiver noise and beam properties, are available from 
the Legacy Archive for Microwave Background Data Analysis 

1 LAMBDA) websi te^. The maps are provided in the HEALPix"^ 
Gorski et all 1999!) format at a resolution of N side = 512 (the num- 
ber of pixels in a HEALPix m ap is given by 12 Nsm /). The data 
processing pipeline specified bv'Ko matsu et alJi2003 ^ is applied to 
produce a single co-added map for analysis. The co-added temper- 
ature at a given position on the sky oj is given by 



T(0J): 






(11) 



where T{(jj) is a CMB temperature map and the r index corresponds 
to the Q-, V- and W-band receivers respectively (indices r = 1,2 
correspond to the K and Ka receiver bands that are excluded from 
the analysis). The noise weights Wrico) are defined by 



Wr{(jS) ■■ 






(12) 



^ http://cmbdata.gsfc.nasa.gov/ 

^ http ://w w w. eso . org/science/healpix/ 
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(a) Mexican hat e = 0.00 (b) Mexican hat e = 0.95 (c) Real Morlet k = (10, 0)^ 

Figure 2. Spherical wavelets (dilation a = 750'; size on sky ^i = 2100', ^2 = 236'). 



where Nrico) specifies the number of observations at each point on 
the sky for each receiver band and ctq^ is the receiver noise disper- 
sion. 

Foreground cleaned sky maps, where the Galactic foreground 
signal (consisting of synchrotron, free-free, and dust emission) has 
been removed, are directly available from the LAMBDA website. 
The Galactic foreground signal is removed by usin g the 3 -band, 
5 -para meter template fitting method described by iBennett et alJ 



y003|). We use these foreground cleaned maps in our analysis. 

An in dependen t foreground analysis of the WMAP data is per- 
formed by Tesmark et al. (2003). The Tegmark et al. map is also 
constructed from a linear summation of observed WMAP maps, 
however the weights used vary over both position on the sky and 
scale. We also perform our analysis on the Tegmark et al. map to 
ensure any detected deviations from Gaussianity are not due to dif- 
ferences in the various foregr ound removal techniques. 

Following the analysis o rfVielva et al J 12003 ) we down-sample 
map resolutions to Nstde = 256, since the very small scales are dom- 
inated by noise (and also to reduce computational requirements). 
The conservative KpO exclusion mask provided by the WMAP 
team (appropriately downsampled to conserve point source exclu- 
sion regions in the coarser resolution) is applied to remove emis- 
sions due to the Galactic plane and known point sources. The final 
preprocessed co-added map (hereafter referred to as the WMAP 
team map, or simply WMAP map) and the map produced by 
iTegmark et alJ i2003h (hereafter referred to as the Tegmark map) 
are illustrated in Figure El 



3.2 Monte Carlo simulations 

Monte Carlo simulations are performed to construct confidence 
bounds on the test statistics used to probe for non-Gaussianity in 
the WMAP 1-year data. 1000 Gaussian CMB realisations are pro- 
duced from the theoretical power spectrum fitted by the WMAP 
team.^ 

To simulate the WMAP observing strategy each Gaussian 
CMB realisation is convolved with the beam transfer function of 
each of the Q-, V- and W-band receivers. White noise of disper- 
sion cr f.(cL>) - (To J ^JNr(oJ) is added to each band. The resultant 
simulated Q-, V- and W-band maps are combined in the same man- 
ner used to construct the co-added map, before down-sampling and 



^ The theoretical power spectrum used is based on a Lambda Cold Dark 
Matter (ACDM) model using a power law for the primordial spectral in- 
dex which best fits the WMAP, Cosmic Background Imager (CBI) and Ar- 
cminute Cosmology Bolometer Array Receiver (ACBAR) CMB data, and 
is also directly available from the LAMBDA website. 



applying the KpO mask, to give a final simulated Gaussian co-added 
map for analysis. 

The same Gaussian simulations are also used for compari- 
son with the Tegmark map. Since the weights used to construct 
the Tegmark map diff'er from those used to construct the WMAP 
team map, one should strictly produce a second set of Gaussian 
simulations following the Tegmark map construction method. The 
weights for the Tegmark map vary as a function of angular scale, 
and unfortunately are not quoted explicitly. Nevertheless, for both 
the WMAP and Tegmark maps, the weights sum to unity and the 
diff'erence in the linear combination of maps used by Teg mark et alJ 
(2003) should not lead to significant changes in the Gaussian con- 
fidence limits as compared with those obtained using the Gaussian 
simulations produced to model the WMAP map. 



3.3 Wavelet analysis 

The CSWT is a linear operation; hence the wavelet coefficients of 
a Gaussian map will also follow a Gaussian distribution. One may 
therefore probe a full sky CMB map for non-Gaussianity simply 
by looking for deviations from Gaussianity in the distribution of 
the spherical wavelet coefficients. 

The analysis consists of first taking the CSWT at a range of 
scales and, for directional wavelets on the sphere, a range of 7 
directions. The scales we consider (and the corresponding effective 
size on the sky for both the Mexican hat and real Morlet wavelets) 
are shown in Table [l] For directional wavelets we consider five 
evenly spaced 7 orientations in the domain [0, tt). 

Those wavelet coefficients distorted by the application of the 
KpO mask are removed, as subsequently described, before test 
statistics are calculated from the wavelet coefficients. An identi- 
cal analysis is performed on each Monte Carlo CMB simulation in 
order to construct significance measures for the test statistics. 



3. 3. 1 Coefficient exclusion masks 

The application of the KpO exclusion mask distorts coefficients cor- 
responding to wavelets that overlap with the mask exclusion region. 
These contaminated wavelet coefficients must be removed from any 
subsequent non-Gaussianity analysis. An extended coefficient ex- 
clusion mask is required to remove all contaminated wavelet coef- 
ficients. 

On small scales masked point sources introduce significant 
distortion in wavelet coefficient maps and should not be neglected. 
On larger scales the masked Galactic plane introduces the most sig- 
nificant distortion, as point source distortions are averaged over 
a large wavelet support. Our construction of an extended coeffi- 
cient mask inherently accounts for the dominant type of distortion 
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(a) Co-added WMAP team map (b) Tegmark map 

Figure 3. Preprocessed WMAP maps considered in the non-Gaussianity analysis. 



Table 1. Wavelet scales considered in the non-Gaussianity analysis. The overall size on the sky ^i for a given scale are the same 
for both the Mexican hat and real Morlet wavelets. The size on the sky of the internal structure of the real Morlet wavelet ^2 is 
also quoted. 



Scale 



1 



10 



11 



12 



Dilation a 


50' 


100' 


150' 


200' 


250' 


300' 


Size on sky ^\ 


UV 


282' 


424' 


565' 


706' 


847' 


Size on sky ^2 


15.1' 


31.4' 


47.1' 


62.8' 


78.5' 


94.2' 



350' 



110' 



400' 


450' 


500' 


550' 


600' 


1130' 


1270' 


1410' 


1550' 


1690' 


126' 


141' 


157' 


173' 


188' 



on a particular scale. Firstly, the CSWT of the original KpO 
mask is taken. Adm issible spherical wavelets have zero mean 
JAntoine et all2002n . hence the only non-zero wavelet coefficients 
are those that are distorted by the mask boundary. These distorted 
coefficients may be easily detected and the coefficient exclusion 
mask extended accordingly. Coefficient exclusion masks are illus- 
trated in FigureElfor the Mexican hat e - 0.00 wavelet for a range 
of scales and in Figure |5] for the real Morlet wavelet for a given 
scale (the scale that a significant non-Gaussianity detection is sub- 
sequently made) and a range of orientations. As scale increases the 
dominant form of distortion may be seen in Figure |3]to shift from 
poin t source to G a lactic plane. 

IVielva et aP i2003h construct an extended coefficient mask 
simply by extending the Galactic plane region of the KpO mask 
by 2.5a (the point source components of the original mask are not 
extended). Several othe r definitions for coefficient e xclusion masks 
are analysed in detail bv lMukheriee & WangH2004 h. none of which 
alter the results of subsequent non-Gaussianity analysis. Although 
it is important to account correctly for the distortions introduced 
by the KpO mask, the results of Gaussianity analysis appear to be 
relatively insensitive to the particular mask chosen. 



3.3.2 Test statistics 

The third (skewness) and fourth (kurtosis) moments about the mean 
are considered to test spherical wavelet coefficients for deviations 
from Gaussianity. These estimators describe the degree of symme- 
try and the degree of peakedness in the underlying distribution re- 
spectively. Skewness is defined by 



and excess kurtosis by 



aa,y) 



1 



^ [Wt{a,y)-li{a,y)f 
^ cr^ (a, 7) 



(13) 



. ^eff r 



Wi{a,y) - id{a,y)f 



o^{a,y) 



- 3. 



(14) 



where ji is the mean and cr the dispersion of the wavelet coeffi- 
cients. The / index ranges over all wavelet coefficients not excluded 
by the coefficient exclusion mask and indexes both a and p com- 
ponents. The number of spherical wavelet coefficients retained in 
the analysis after the application of the coefficient exclusion mask 
is given by A^efF- 

Skewness and excess kurtosis for a Gaussian distribution are 
both zero. We look for deviations from zero in these test statistics 
to indicate the existence of non-Gaussianity in the distribution of 
spherical wavelet coefficients, and hence also in the corresponding 
CMB map. 



4 RESULTS 

To probe for non-Gaussianity in the WMAP 1-year data, the anal- 
ysis procedure described in Section |3l is performed on both the 
WMAP team and Tegmark maps. The three spherical wavelets il- 
lustrated in Figure |3 are considered, namely the symmetric Mex- 
ican hat e - 0.00 wavelet, the elliptical Mexican hat e - 0.95 
wavelet and the real Morlet k - (10, 0)^ wavelet. The Mexican hat 
6 = 0.00 case has previously been analysed by Vielva et al. (2003) 
(although some scales considered diff'er), thereby providing a con- 
sistency check for the analysis. 



4.1 Wavelet coefficient statistics 

For a given wavelet, the skewness and kurtosis of wavelet 
coefficients is calculated for each scale and orientation. These 
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(a) ai = 50' 



(b) «2 = 100' 



(c) as = 150' 




(d) (34 = 200' 



(e) as = 250' 



(f) ae = 300' 




(g) aj = 350' 



(h) as = 400' 



(i) fl9 = 450' 




(j) flio = 500' (k) an = 550' (1) an = 600' 

Figure 4. Symmetric Mexican hat e = 0.00 wavelet coefficient exclusion masks for each scale. 




(a) 7 = 0° 



(b) r = 72° 



(c) 7 = 144° 




(d) 7 = 216° (e) 7 = 288° 

Figure 5. Real Morlet wavelet coefficient exclusion masks at scale a\\ = 550' for each orientation. 
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statistics are displayed in Figure |6l with confidence intervals 
constructed from the Monte Carlo simulations also shown. For 
directional wavelets, only the orientations corresponding to the 
maximum deviations from Gaussianity are shown. 

Our coefficient e x clusio n mask diff'ers slightly from that ap- 
plied by IVielva et all i2003h . thus for comparison purposes we 
also perform the Mexican hat e - 0. 00 analysis withou t apply- 
ing any extended coefficient mask, as Vielva et al. ( 2003) also do 
initially. These results, although not shown, correspond identically. 
By applying diff'erent coefficient masks the shape of the plots diff'er 
slightly, nevertheless the findings drawn remain the same. Devia- 
tions from Gaussianity are detected in the kurtosis outside of the 
99% confidence region constructed from Monte-Carlo simulations, 
on scales a^ - 250' and ^6 = 300'. Furthermore, a deviation out- 
side the 9 9% confidence regio n is detected in the skewness at scale 
^2 = 100'. I Vielva et alJf2003l) measure a similar skewness value at 
this scale, although this lies directly on the boundary of their 99% 
confidence region. 

Deviations from Gaussianity are also detected in both skew- 
ness and kurtosis using the elliptical Mexican hat e - 0.95 wavelet. 
In each case the observed deviations occur on a slightly larger 
scale than those found using the symmetric Mexican hat e - 0.00 
wavelet. This behaviour also appears typical for simulated Gaus- 
sian map realisations. Adjacent orientations exhibit similar results, 
although not at such large confidence levels (but still outside of the 
99% confidence level). 

An extremely significant deviation from Gaussianity is ob- 
served in the skewness of the real Morlet wavelet coefficients at 
scale a\i - 550' and orientation y = 72°. The kurtosis measure- 
ment on the same scale and orientation also lies outside of the 99% 
confidence region. 



4.2 Statistical significance of detections 

We now consider in more detail, the most significant deviation from 
Gaussianity obtained in each of the panels in Figure |6l In particu- 
lar, we examine the distribution of each statistic, obtained from the 
Gaussian Monte Carlo simulations, and also perform a x^ test for 
each statistic. Significance measures of the non-Gaussianity detec- 
tions may then be constructed from each test. 

FigureEl shows histograms constructed from the Monte Carlo 
simulations for those test statistics corresponding to the most sig- 
nificant deviations from Gaussianity. The measured statistic for 
both the WMAP team and Tegmark maps is also shown on each 
plot, with the number of standard deviations these observations de- 
viate from the mean. In particular, we note the large deviations 
shown in panel (c), corresponding to 5.61 and 6.42 standard de- 
viations for the real Morlet wavelet analysis of the WMAP and 
Tegmark maps respectively. 

Having determined separately the confidence level of the 
largest non-Gaussianity detection in each panel of Figure|6l we now 
consider the statistical significance of our results for each wavelet 
as a whole. Treating each wavelet separately, we search through the 
Gaussian simulations to determine the number of maps that have 
an equivalent or greater deviation in any of the test statistics calcu- 
lated from that map using the given wavelet. That is, if any skew- 
ness or kurtosis statistic^ calculated from the Gaussian map - on 



^ Although we recognise the distinction between skewness and kurtosis, 
there is no reason to partition the set of test statistics into skewness and 
kurtosis subsets. The full set of test statistics must be considered. 



Table 2. Deviation and significance levels of spherical wavelet coefficient 
statistics calculated from the WMAP map (similar results are obtained using 
the Tegmark map). Standard deviations and significant levels are calculated 
from 1000 Monte Carlo Gaussian simulations. The table variables are de- 
fined as follows: the number of standard deviations the observation deviates 
from the mean is given by Ng-', the number of simulated Gaussian maps that 
exhibit an equivalent or greater deviation in any test statistics calculated us- 
ing the given wavelet is given by A^^dev; the corresponding significance level 
of the non-Gaussianity detection is given by 6. Only those scales and ori- 
entations corresponding to the most significant deviations from Gaussianity 
are shown for each wavelet. 



Skewness 
(«2 = 100') 



Kurtosis 
{ae = 300') 



A^dev 

s 



-3.38 
28 maps 

97.2% 



3.12 

47 maps 

95.3% 



(a) Mexican hat e = 0.00 



Skewness Kurtosis 

(fl3 = 150'; 7 = 72°) (^lo = 500'; y = 108°) 



A^dev 

5 



-4.10 
39 maps 
96.1% 



3.01 

199 maps 

80.1% 



(b) Mexican hat e = 0.95 



Skewness Kurtosis 

(an = 550'; y = 72°) (an = 550'; y = 72°) 



A^dev 

5 



-5.61 
17 maps 
98.3% 



2.66 
642 maps 

35.8% 



(c) Real Morlet A: = (10,0)^ 



any scale or orientation - deviates more than the maximum devi- 
ation observed in the WMAP data for that wavelet, then the map 
is flagged as exhibiting a more significant deviation. This is an ex- 
tremely conservative means of constructing significance levels for 
the observed test statistics. Significance levels corresponding to the 
detections considered in Figure are calculated and displayed in 
Table |2| Notice that although several individual test statistics fall 
outside of the 99% confidence region, the true significance level 
of the detection when all statistics are taken into account is con- 
siderably lower. Using our conservative test, the significance of the 
non-Gaussianity detection made by Vielva et al. (2003), previously 
quoted at greater than 99.9% significance, drops to a significance 
level of 95.3%. Of particular interest is the non-Gaussian detec- 
tion in the skewness of real Morlet wavelet coefficients on scale 
an = 550' and orientation y = 72°. This statistic deviates from the 
mean by 5.61 standard deviations for the WMAP map and by 6.42 
standard deviations for the Tegmark map. The detection is made at 
an overall significance level of 98.3%. 

The preceding analysis is based on the marginal distributions 
of individual statistics and makes a posterior selection of the critical 
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(d) Kurtosis - Mexican hat e = 0.00 




(e) Kurtosis 
7 = 108° 
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(f) Kurtosis - real Morlet A: = (10,0)^; 

7 = 72° 



Figure 6. Spherical wavelet coefficient statistics for each wavelet. Confidence regions obtained from 1000 Monte Carlo simulations are shown for 68% 
(red), 95% (orange) and 99% (yellow) levels, as is the mean (solid white line). Only the orientations corresponding to the most significant deviations from 
Gaussianity are shown for the Mexican hat e = 0.95 and real Morlet wavelet cases. 



confidence limit from the most discrepant values obtained from 
the data. A x^ test provides an alternative analysis and method of 
constructing significance measures. This test instead considers the 
set of test statistics for each wavelet as a whole and hence is based 
on their joint distribution. The posterior statistic selection problem 
is thus eliminated, however including a large number of less use- 
ful test statistics has a pronounced efifect on down-weighting the 
overall significance of the test. The;^^ statistic is given by 



AW AW 



X 



J^Yj(r,-Td{C-%{Tj-Tj), 



(15) 



i=i j=i 



where r^ gives each test statistic. For Gaussian distributed test 
statistics this should satisfy a x^ distribution. Although our test 
statistics are not Gaussian distributed, one may still use the x^ test 
if one is willing to estimate significance levels using Monte Carlo 
simulations. The test statistics include both skewness and kurtosis 
statistics for each scale and orientation, hence there are 24 statis- 
tics for the symmetric Mexican hat analysis and 120 for each of the 
directional wavelet analyses. The mean value for each test statis- 
tic T and the covariance matrix C is calculated from the Gaussian 
simulations, x^ values are calculated for the WMAP map, and also 
for all simulated Gaussian realisations. The previously described 
approach for constructing significance levels is applied to the x^ 
statistics. Figure [Sl summarises the results obtained. All spherical 
wavelet analyses flag deviations from Gaussianity of very high sig- 
nificance when all test statistics are incorporated in this manner. 
In particular, the detection made using the symmetric Mexican hat 



wavelet occurs at the 99.9% significance level and that made with 
the real Morlet wavelet occurs at the 99.3% significance level. In 
this case the superiority of the symmetric Mexican hat wavelet over 
the real Morlet wavelet arises since the real Morlet wavelet analy- 
sis contains a number of additional less useful statistics (due to the 
additional orientations), that dilute the overall results. 

Since the first analysis is based on marginal distributions, as 
opposed to the joint distribution of statistics probed by the;^^ test, 
the former is more conservative. Thus we quote the overall signif- 
icance of all detections of non-Gaussianity at the level calculated 
by the first method (which, in all cases, is the lower of the values 
calculated by the two methods). 



4.3 Localised deviations from Gaussianity 

Wavelet analysis inherently aflfords the spatial localisation of 
interesting signal characteristics. The most pronounced deviations 
from Gaussianity in the WMAP 1-year data may therefore be loca- 
lised on the sky. In addition, directional wavelets also allow signal 
components to be localised in orientation. 

The wavelet coefficients corresponding to the most signifi- 
cant non-Gaussian detections for each wavelet are displayed in 
Figure 121 accompanied by corresponding thresholded maps to lo- 
calise the most pronounced deviations from Gaussianity. The re- 
gions displayed in Figure |9l (b) that are detected from the kurto- 
sis Mexican hat 6 = 0.00 analysis are in close accordance with 
those regions found bv lVielva et alJ i2003l ). Additional similarities 
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Figure 7. Histograms of spherical wavelet coefficient statistic obtained from 1000 Monte Carlo simulations. The mean is shown by the dashed vertical line. 
The observed statistics for the WMAP and Tegmark maps are shown by the blue and green lines respectively. The number of standard deviations these 
observations deviate from the mean is also displayed on each plot. Only those scales and orientations corresponding to the most significant deviations from 
Gaussianity are shown for each wavelet. 



appear to exist between the regions detected from different thresh- 
olded wavelet coefficient maps, as apparent in Figure |9l To quan- 
tify these similarities, the cross-correlation of all combinations of 
thresholded coefficient maps is computed (the cross-correlation is 
normalised to lie in the range [-1,1], where unity indicates a fully 
correlated map). Table [S] shows the normalised cross-correlation 
values obtained. Deviation regions shown in Figure |9l(b) and Fig- 
ure 121 (d), detected by the Mexican hat e = 0.00 and e = 0.95 
wavelets respectively, are highly correlated. Interestingly, these co- 
efficient maps are both flagged by excess kurtosis measures. Fur- 
thermore, the regions shown in Figure |9l (a) and Figure |9l(c), de- 
tected by the Mexican hat € = 0.00 and € = 0.95 wavelets respec- 
tively, are moderately correlated. These coefficient maps are both 
flagged by excess skewness measures. No other combinations of 
thresholded coefficient maps exhibit any significant similarities. In 
particular, the deviation regions detected by the skewness real Mor- 
let analysis (Figure |9](e)) do not correlate with any of the regions 
found using the Mexican hat wavelets. This is expected since a dif- 
ferent wavelet that probes different structure is applied. The cross- 
correlation relationships exhibited here between detected deviation 
regions for the WMAP map, also appear typical of simulated Gaus- 
sian maps. 

To investigate the impact of these localised regions on the ini- 
tial non-Gaussianity detection, the corresponding coefficients are 
removed from the calculation of skewness and kurtosis test statis- 
tics. The non-Gaussian detections are substantially reduced for all 



of the six most significant test statistics considered in Figure^l For 
the statistics considered in Figure 0(c), (d) and (e) the detection of 
non-Gaussianity is completely eliminated. For the remaining cases 
considered in FigureEl(a), (b) and (f) non-Gaussian detections are 
reduced in significance and lie between the 95% and 99% confi- 
dence levels. 

Thus, the localised deviation regions identified do indeed ap- 
pear to be the source of detected non-Gaussianity. Moreover, those 
detected regions shown in Figure |9l (a), (c) and (e) appear to 
introduce skewness into the WMAP map, whereas those detected 
regi ons shown i n Figure El (b) and (d) appear to introduce kurtosis. 

ICruz et alJ i2004h . in a continuation of the work of 
IVielva et aP |20pj), consider localised regions in more detail and 
find the large cold spot at (b = -57°, / = 209°) to be the main 
source of non-Gaussianity detected by the symmetric Mexican hat 
wavelet analysis. A similar analysis may be performed using the 
wavelets we consider and the associated localised regions, although 
we leave this for a further work. 



4.4 Preliminary noise analysis 

Naturally, one may wish to consider possible sources of the non- 
Gaussianity detected. We briefly consider the deviation regions de- 
tected to see if they correspond to regions on the sky that have 
higher noise dispersion than typical. We leave the analysis of resid- 
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Figure 9. Spherical wavelet coefficient maps (left) and thresholded maps (right). To locaHse most likely deviations from Gaussianity on the sky, the coefficient 
maps exhibiting strong non-Gaussianity are thresholded so that only those coefficients above 3cr (in absolute value) are shown. Due to the apparent similarity 
of the WMAP team and Tegmark maps, only coefficients for the analysis of the WMAP map are shown above. 
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Figure 8. Histograms of normahsed ;t'^ test statistics obtained from 1000 
Monte Carlo simulations. The normalised x'^ value obtained from the 
WMAP map is indicated by the blue vertical line (similar results are ob- 
tained using the Tegmark map). The number of simulated maps that exhibit 
a greater or equivalent x'^ value than the WMAP map is quoted (Ne), ac- 
companied by the corresponding significance level (SL). 



Table 3. Normalised cross-correlation of thresholded spherical wavelet co- 
efficient maps indicating the similarity between the localised most likely 
deviations from Gaussianity flagged by the most significant skewness and 
kurtosis observations for each wavelet. Notice that the regions detected 
from the skewness flagged maps of the symmetric and elliptical Mexican 
hat wavelets are moderately correlated, while the regions detected from the 
kurtosis flagged maps of the symmetric and elliptical Mexican hat wavelets 
are strongly correlated. The regions flagged by the real Morlet wavelet anal- 
ysis are not correlated with any of the other regions detected by a Mexican 
hat wavelet analysis, as expected since a different wavelet that probes dif- 
ferent structure is applied. (Note that the lettered key corresponds to the 
thresholded coefficient maps contained in the panels of Figure|9l) 
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ual foregrounds or further systematics, or whether the features de- 
tected do indeed exist in the CMB, for a further work. 

The noise dispersion map for each WMAP band is combined, 
according to 



o-(aj) : 






(16) 



to produce a single noise dispersion sky map for the WMAP map, 
and equivalently for the Gaussian realised maps. A histogram of 
this map, and a KpO masked version of the map, is illustrated in 
Figure [TOI to investigate the noise dispersion distribution for the 
WMAP observing strategy. Also plotted is the mean noise disper- 
sion level in the detected deviation regions for each thresholded 
coefficient map of Figure|9l All mean noise dispersion levels of de- 
viation regions lie within the central region of the full noise distri- 
bution. Furthermore, full noise dispersion histograms for the devia- 
tion regions were also produced to ensure outliers did not exist. No 
outliers were observed in any deviation regions. These additional 
five plots are not shown to avoid clutter and since no pertinent ad- 
ditional findings may be drawn from them. It is therefore apparent 
that the deviation regions detected do not correspond to regions 
with greater noise dispersion than typical. 



5 CONCLUSIONS 

A directional spherical wavelet analysis, facilitated by our fast 
CSWT, has been applied to the WMAP 1-year data to probe for de- 
viations from Gaussianity. Directional spherical wavelets allow one 
to probe orientated structure inherent in the data. Non-Gaussianity 
has been detected by a number of test statistics for a range of 
wavelets. 

W e have reproduced the results obtained by IVielva et all 
i2003h using the symmetric Mexican hat 6 = 0.00 wavelet, thereby 
confirming their findings, whilst also providing a consistency check 
for our analysis. Deviations in the skewness and kurtosis of wavelet 
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Figure 10. WMAP sky noise dispersion histogram, with the mean noise dis- 
persion level obtained in detected deviation regions also shown. The solid 
histogram shown corresponds to the full sky noise dispersion map, whereas 
the dotted histogram corresponds to the KpO masked noise dispersion map. 
The dashed vertical lines indicate the mean noise dispersion level in the 
detected deviation regions illustrated in Figure |9| None of the detected de- 
viation regions correspond to noise of higher dispersion than typical. The 
lettered key corresponds to the thresholded coefficient maps contained in 
the panels of Figure |9l (a) Mexican hat e = 0.00, a2 = 100'; (b) Mex- 
ican hat € = 0.00, ae = 300'; (c) Mexican hat e = 0.95, as = 150', 
y = 12°; (d) Mexican hat e = 0.95, aio = 500', y = 108°; (e) Real Morlet 
k = im,Of,an =550', 7 = 72°. 



coefficients on scale ^2 = 100' and a^ - 300' were detected, al- 
though using our more conservative test we make these detections 
at the 97.2% and 95.3% overall significance levels respectively 
(lower than the 99.9% significance level quoted by IVielva et alJ 
J2QQ3h for their kurtosis detection). 

Similar detections of non-Gaussianity were made using the 
elliptical Mexican hat e = 0.95 wavelet, although on slightly larger 
scales. In particular, a deviation from Gaussianity was detected in 
the skewness of the Mexican hat € = 0.95 wavelet coefficients on 
scale as = 150' and orientation y = 72° at the 96.1% significance 
level. Although a detection was observed in the kurtosis outside 
of the 99% confidence region on scale aio = 500' and orientation 
7 = 108°, the full statistical analysis of Monte Carlo simulations 
gave a significance of only 80. 1 % for this detection. 

The most interesting result, however, is the deviation from 
Gaussianity observed in the real Morlet wavelet skewness measure- 
ment on scale an = 550' and orientation y = 72°. This wavelet 
scale corresponds to an eff'ective size on the sky of ~ 26° (~ 3° for 
the internal structure of the real Morlet wavelet), or equivalently a 
spherical harmonic scale of ^ ~ 7 (^ ~ 63). The detection deviates 
from the mean of 1000 Gaussian Monte Carlo simulations by 5.61 
standard deviations for the WMAP map and by 6.42 standard devi- 
ations for the Tegmark map. Only 17 of 1000 Gaussian simulated 
maps exhibited a deviation this large in any real Morlet test statistic, 
hence the detection is conservatively made at 98.3% significance. 

Significance levels were also calculated from;^;^ tests for each 
spherical wavelet. This approach avoids the posterior selection of 
particular statistics, but rather considers the set of test statistics in 
aggregate. By considering the joint distribution of test statistics in 
this manner the analysis results may be diluted by including a large 
number of less powerful test statistics. Deviations from Gaussianity 
at significance levels of 99.9% and 99.3% were found using the 
symmetric Mexican hat and real Morlet wavelet respectively. In this 



case the directional real Morlet analysis is more severely affected 
by a larger number of less useful test statistics, nevertheless both 
deviations from Gaussianity are made at very high significance. We 
quote the overall significance of our findings, however, at the lower 
significance levels found using the previous most conservative test. 

Deviations from Gaussianity corresponding to the most sig- 
nificant detections for each wavelet were localised on the sky. 
By removing the coefficients corresponding to these regions from 
the initial analysis, all significant non-Gaussianity detections were 
eliminated. These localised regions therefore appear to be the 
source of detected non-Gaussianity. Moreover, those regions that 
introduce skewness in the WMAP map may be localised, as may 
those regions that introduce kurtosis. Preliminary noise analysis in- 
dicates that these detected deviation regions do not correspond to 
regions that have higher noise dispersion than typical. Further anal- 
ysis is required, however, to ascertain whether these regions corre- 
spond to the localised introduction of secondary non-Gaussianity 
or systematics, or whether in fact the non-Gaussianity detected in 
the WMAP 1-year data is due to intrinsic primordial fluctuations in 
the CMB. 

An interesting first step in deducing whether the non-Gaussian 
signal discovered is of cosmological origin would be to repeat our 
analysis on the 4-year COBE-DMR data. Although it has been 
shown that these data contain some systematic effects that lead 
to non-Gaussianity (Masueiio & Medeiros 2004), it is likely that 
these systematics are not shared by WMAP. Provided the most sig- 
nificant detection in the WMAP data is predominantly due to the 
global structure of the real Morlet wavelet at an angular scale of 
~ 26°, the angular resolution of the COBE-DMR data should be 
sufficient to observe it if it is astrophysical in origin. Clearly, it will 
also be of great interest to investigate whether the non-Gaussianity 
detections reported here are still present in the 2-year WMAP data. 
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APPENDIX A: A FAST DIRECTIONAL CSWT 

The CSWT at a particular scale is essentially a spherical convolu- 
tion; hence it is possib le to apply the fas t spher ical convolution 
algorithm proposed by IWandelt & GorskJ i200ll) to evaluate the 
wavelet transform. The harmonic representation of the CSWT is 
first presented, followed by a discretisation and fast implementa- 
tion. Without loss of generality we consider a single dilation only. 



Al Harmonic formulation 

There does not exist any finite point set on the sphere that is invari- 
ant under rotations (due to geometrical properties of the sphere). 



hence it is more natural, and in fact more accurate for numerical 
purposes, to recast the CSWT defined by (|3j in harmonic space. 

Both the wavelet and signal are represented in terms of a 
spherical harmonic expansion, defined for an arbitrary function 



/(^) = Z Z /^-^^-(^) ' 



(Al) 



where the spherical harmonic coefficients are given by the usual 
projection of the signal onto each spherical harmonic basis function 



h 



J52 



(A2) 



In practice one requires that at least one of the functions, usually the 
wavelet, has a finite band limit so that negligible power is present 
in those coefficients above a certain ^^ax- For all practical purposes, 
the outer summation of (lAU may then be truncated to ^max- 

Substituting the spherical harmonic expansions of the wavelet 
and signal into the CSWT of (|3j and noting the orthogonality of the 
spherical harmonics, yields the harmonic representation 



W(a,(i,y) = J] J] J] [Di„,{a^l3,y)^en,'] 



Sim 



(A3) 



e=Q m=-e m'-- 



The additional summation and Z)^^^^, Wigner rotation matrices that 
are introduced characterise the rotation of a spherical harmonic, 
noting that a rotated spherical harmonic may si mply be represent ed 
by a sum of rotated harmonics of the same i by fcnui et alll996h 

I 



The Wigner rotation matrices may be decomposed as 



DLA»'P^y)- 



'dL'We- 



(A4) 



(A5) 



where the real polar J-matrix is defined by, for example, 
iBrink & Satchlei| il993h . The relationship shown in (IA5t is ex- 
ploited by factoring the rotation ^a,fi,y into two separate rotations, 
both of which only contain a constant ±7t/2 polar rotation: 



'^a,l3,y - % 



■a-n/2,-n/2,/3 



%>. 



n/2, y+n/2 • 



(A6) 



By factoring the rotation in this manner and applying the decom- 
position described by (IA5L (IA3t can be written as 



^max £ i min(mmax/) 

W(a,li,y) = Y,Y. Y. E d.i,„,(nl2)di,„,,{7'l2) 

£=0 m=-£ m'=-t m"=-min(mmax/) 

where the symmetry relationship S^^X-p) - ^t'm^^ ^^^ ^^^^ ^P" 
plied. In many cases it is likely that the wavelet will have mini- 
mal azimuthal structure compared to the signal under analysis, in 
which case it may also have a lower effective azimuthal band limit 

'^max "^ ^max' 

The harmonic formulation presented replaces the continuous 
integral of (|3} by finite summations, although evaluating these sum- 
mations directly would be no more efficient that approximating the 
initial integral using simple quadrature. Rotations are elegantly rep- 
resented in harmonic space, however, and the approximation and 
interpolation required in any real space discretisation is avoided. 
Moreover, (IA7t is represented in such a way that the presence of 
complex exponentials may be exploited such that fast Fourier trans- 
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forms (FFTs) may be applied to evaluate rapidly the three summa- 
tions simultaneously. 



A2 Fast implementation 

Azimuthal rotations may be applied with far less computational ex- 
pense than polar rotations since they appear within complex expo- 
nentials in \A1\ . Although the J-matrices can be evaluated reason- 
ab ly quickly an d reliably using recursion formulae (e.g. those given 
bv lRisbdll996l) . the basis for the fast implementation is to avoid 
these polar rotations as much as possible and use FFTs to evalu- 
ate rapidly all of the azimuthal rotations simultaneously. This is the 
motivation for factoring the rotation by (IA6t so that all Euler angles 
occur as azimuthal rotations. 

The discretisation of each Euler angle may in general be arbi- 
trary. However, to exploit standard FFT routines uniform sampling 
is adopted. The uniformly sampled spherical wavelet coefficient 
samples are defined by'' 

/ Inricy InriR Inriy \ 
^...... = ^(^.^.^)- (A8) 

Discretising (IA7i in this manner and performing a little algebra we 
may recast it in a form amenable to fast Fourier techniques: 

W_ -i2n(na£max/Na+nBemax/NB+nymmax/Ny) 
na,nB,ny — t- 

Na-l ^yS-1 Ny-l 
X V V \^ ^. .^ .^^ ^'^^UnalNa+fnplNp+f'nylNy) ^ /^^x 

j=0 /=0 /'=0 

where the second line is simply the unnormalised 3D inverse dis- 
crete Fourier transform (DFT) of 

f _ i(m"-m)7T/2 

'-m+imaxM +€max,in" +mmax ~ ^ 

Anax 

X X «'».'m('^/2) €'„," (.^1^) ^f"" * ^fm - (AlO) 

^=max(|m|,|OT'|,|m"|) 

where the shifted indices show the conversion between the 
harmonic and Fourier conventions. The number of samples 
for each Euler angle is A^^^ = 2^max + 1, ^p - 2^max + 1 and 
Ny - 2mmax + 1, enforced by uniform sampling and the standard 
Fourier relationship. 

The CSWT may be performed very rapidly in spherical har- 
monic space by constructing the ^matrix of JAlOt from spher- 
ical harmonic coefficients and precomputed J-matrices, followed 
by the application of a FFT to evaluate rapidly all three Euler an- 
gles of the discretised CSWT simultaneously, before applying a fi- 
nal modulating complex exponential factor. Memory and compu- 
tational requirements may be reduced by a further factor of two 
for real signals by exploiting the conjugate symmetry relationship 

t-m-m'-m" - ^m,m',m"- 

The computational cost of the fast CSWT is dominated by 
the calculation of the ^matrix, which is of order Oii^^^m^^y^) . 
A direct quadrature approximation of the CSWT integral is of 
order O{Na^N0^Ny). The harmonic and real space size parame- 
ters are of the same order, that is (9(^inax) ~ 0(Na) ~ C)(A/^) and 
0(mjnax) ~ 0(Ny), hence the fast algorithm provides a saving of 
<^(^max)- We give a more detailed comparison of the complexity 
of various CSWT ini plementations and typical execution times in 
iMcEwen et"allj2004 . 



^ Whilst a and y both cover the range to In, evaluating j3 over the same 
range is redundant, covering the S0(3) manifold exactly twice. Nonethe- 
less, the use of the fast FFT-based algorithm requires this range. 



